
#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern int number_time_steps;

#define max_value 1000
#define dsource(ix,iy) gddd(dd,(ix),(iy))
#define dsource1(incx,incy) dsource(ix1+(incx),iy1+(incy))
#define ddest(ix,iy) gddd(dd1,(ix),(iy))
#define dtest(ix,iy) gddd(d4,(ix),(iy))

#define ckchange(p1,p2) \
  ( ((p1<=0.0)&&(p2>0.0))||((p1>0.0)&&(p2<=0.0)) )

#define cknode(ixinc,iyinc) \
  (ckchange(dsource(ix,iy),dsource(ix+ixinc,iy+iyinc)) )

logical transition(ix,iy,dd)
int ix,iy;
double dd[imax][jmax];
{
  if ( (cknode(0,1)) || (cknode(1,0)) || (cknode(-1,0)) || (cknode(0,-1)) )
   return TRUE;
  else
   return FALSE;
}   /* transition */


#define cknode1(ixinc,iyinc) \
 ( ((ddest(ix,iy)<=-max_value)&&(dtest(ix+ixinc,iy+iyinc)>-max_value))|| \
   ((ddest(ix,iy)>=max_value)&&(dtest(ix+ixinc,iy+iyinc)<max_value)) \
 )

logical transition1(ix,iy,dd1)
int ix,iy;
double dd1[imax][jmax];
{
  if ( (cknode1(0,1)) || (cknode1(1,0)) || (cknode1(-1,0)) || (cknode1(0,-1)) )
   return TRUE;
  else
   return FALSE;
}   /* transition1 */

void eval_distance(i1,j1,p1,i2,j2,p2,i3,j3,p3,ix,iy,ix1,iy1,out_dist)
int i1,j1,i2,j2,i3,j3,ix,iy,ix1,iy1;
double p1,p2,p3,*out_dist;
{
double factor1,factor2,dist1,dist2,dist3,
  x1,y1,x2,y2,x3,y3,xn,yn,slope,slope1,b,b1;

  if ( (ckchange(p1,p2))&&(ckchange(p2,p3)) )
   {
    if (fabs(p1)+fabs(p2)>0.0)
     factor1=fabs(p1)/(fabs(p1)+fabs(p2));
    else
     factor1=0.5;

    if (fabs(p2)+fabs(p3)>0.0)
     factor2=fabs(p2)/(fabs(p2)+fabs(p3));
    else
     factor2=0.5;

    x1=h*(ix1+factor1*(i2-i1)+i1);
    y1=h*(iy1+factor1*(j2-j1)+j1);
    x2=h*(ix1+factor2*(i3-i2)+i2);
    y2=h*(iy1+factor2*(j3-j2)+j2);
    x3=h*ix;y3=h*iy;

    if (x2==x1)
     {
      xn=x1;
      yn=y3;
     }
    else if (y2==y1)
     {
      xn=x3;
      yn=y1;
     }
    else
     {
      slope=(y2-y1)/(x2-x1);
      b=y1-slope*x1;
      slope1=-1/slope;
      b1=y3-slope1*x3;
      xn=(b-b1)/(slope1-slope);
      yn=slope*xn+b;
     }

    dist1=sqrt(sqr(x1-x3)+sqr(y1-y3));
    dist2=sqrt(sqr(x2-x3)+sqr(y2-y3));
    if ((xn>=min(x1,x2))&&(xn<=max(x1,x2)) )
     dist3=sqrt(sqr(xn-x3)+sqr(yn-y3));
    else
     dist3=max_value;
    dist1=min(dist1,dist2);
    *out_dist=min(dist1,dist3);
   }   /* line went through the triangle */
  else
   *out_dist=max_value;
}  /* eval_distance */

double set_distance(ix,iy,ix1,iy1,dd)
int ix,iy,ix1,iy1;
double dd[imax][jmax];
{
double current_distance,out_dist;

  current_distance=max_value;

#define check_line(i1,j1,i2,j2,i3,j3) \
  eval_distance(i1,j1,dsource1(i1,j1),i2,j2,dsource1(i2,j2), \
     i3,j3,dsource1(i3,j3),ix,iy,ix1,iy1,&out_dist);  \
  if (out_dist<current_distance) current_distance=out_dist;

  /* check for curves in upper right triangle */
  check_line(0,0,1,0,0,1);
  check_line(1,0,0,1,0,0);
  check_line(0,1,0,0,1,0);

  /* check for curves in upper left triangle */
  check_line(0,0,0,1,-1,0);
  check_line(0,1,-1,0,0,0);
  check_line(-1,0,0,0,0,1);

  /* check for curves in lower left triangle */
  check_line(0,0,-1,0,0,-1);
  check_line(-1,0,0,-1,0,0);
  check_line(0,-1,0,0,-1,0);

  /* check for curves in lower right triangle */
  check_line(0,0,0,-1,1,0);
  check_line(0,-1,1,0,0,0);
  check_line(1,0,0,0,0,-1);

  if (current_distance>=max_value)
   printf("distance error-could not find intruding line\n");

  if (dd[ix][iy]<=0.0)
   current_distance=-current_distance;
/*
printf("sdist ix,iy,ix1,iy1,d,dcmp%5d%5d%5d%5d%10.5f%10.5f\n",ix,iy,ix1,iy1,
  current_distance,d1[ix][iy]);
*/
  return current_distance;
}   /* set_distance */

void brute_reset_distance(dd,dd1,avg_error,test_error)
double dd[imax][jmax],dd1[imax][jmax],*avg_error,*test_error;
{
int ix,iy,ix1,iy1,thick_max,thick_count,search_cnt,sum_count,
    ixbase,iybase;
double max_set,min_set,start_d,sign_d,closest,sum_error,dx,dy,
  dist_test,sum_error1;

  for_all_points
   if (point_ok)
    {
     if (transition(ix,iy,dd))
      dd1[ix][iy]=set_distance(ix,iy,ix,iy,dd);
     else if (dd[ix][iy]>0)
      dd1[ix][iy]=max_value;
     else
      dd1[ix][iy]=-max_value;
     d4[ix][iy]=dd1[ix][iy];
    }

  max_set=min_set=0.0;
  thick_max=(int) ((density_spread)+5);
  for (thick_count=1;thick_count<=thick_max;thick_count++)
   {
    for_all_points
     if (point_ok)
      if (transition1(ix,iy,dd1))
        {
         start_d=dd1[ix][iy];
         sign_d=sign(start_d);
         closest=sign_d*max_value;
         ix1=ix;iy1=iy;
         ixbase=ix;iybase=iy;
         search_cnt=1;

#define test_closest(incx,incy,incx1,incy1) \
  if (sign_d*ddest(ixbase+incx,iybase+incy)<sign_d*closest)  \
   { closest=ddest(ixbase+incx,iybase+incy); \
     ix1=ixbase+incx1;iy1=iybase+incy1; }
  
/* 
         printf("d,ix,iy%10.5f%5d%5d\n",start_d,ix,iy); 
*/
         while ((search_cnt<=thick_max+10)&&(sign_d*closest>0.0))
          {
           test_closest(2,0,1,0)
           test_closest(1,1,1,0)
           test_closest(1,-1,1,0)
    
           test_closest(-2,0,-1,0)
           test_closest(-1,1,-1,0)
           test_closest(-1,-1,-1,0)

           test_closest(0,2,0,1)
           test_closest(1,1,0,1)
           test_closest(-1,1,0,1)

           test_closest(0,-2,0,-1)
           test_closest(1,-1,0,-1)
           test_closest(-1,-1,0,-1)

           ixbase=ix1;
           iybase=iy1;
/*
printf("scound,closest,ix,iy,ix1,iy1%5d%10.5f%5d%5d%5d%5d\n",
 search_cnt,closest,ix,iy,ix1,iy1);
*/
           if (sign_d*closest>=max_value)
            printf("error-no closest found\n");

           search_cnt++;
          }    /* looping to the interface */

         if (search_cnt>thick_max+11)
          printf("error-could not get to the interface\n");
         else
          {
           dd1[ix][iy]=start_d=set_distance(ix,iy,ix1,iy1,dd);

           if ((start_d>max_set)&&(start_d<max_value))
            max_set=start_d;
           if ((start_d<min_set)&&(start_d>-max_value))
            min_set=start_d;      
          }
        }  /* un-init point found with init neighbors */

    for_all_points    /* mark new points init during current sweep */
     if (point_ok)
      d4[ix][iy]=dd1[ix][iy];
   }   /* looping thick_count */

  sum_error=0.0;
  sum_error1=0.0;
  sum_count=0;

  for_all_points
   if (point_ok)
    if (dd1[ix][iy]>=max_value)
     dd1[ix][iy]=max_set+h/2.0;
    else if (dd1[ix][iy]<=-max_value)
     dd1[ix][iy]=min_set-h/2.0; 
    else if (fabs(dd1[ix][iy])<(density_spread+1.0)*h)
     {
      dx=(dsource(ix+1,iy)-dsource(ix-1,iy))/(2*h);
      dy=(dsource(ix,iy+1)-dsource(ix,iy-1))/(2*h);
      dist_test=sqrt(sqr(dx)+sqr(dy));
      if (dist_test<almost_same)
       dist_test=1/almost_same;
      else if (dist_test<1.0)
       dist_test=1/dist_test;

      sum_error+=fabs(dist_test-1.0);

      dx=(ddest(ix+1,iy)-ddest(ix-1,iy))/(2*h);
      dy=(ddest(ix,iy+1)-ddest(ix,iy-1))/(2*h);
      dist_test=sqrt(sqr(dx)+sqr(dy));
      if (dist_test<almost_same)
       dist_test=1/almost_same;
      else if (dist_test<1.0)
       dist_test=1/dist_test;

      sum_error1 += fabs(dist_test-1.0);

      sum_count++;
     }
  
  if (sum_count>0)
   {
    *test_error=sum_error1/sum_count;
    *avg_error=sum_error/sum_count;
   }
  else
   {
    *test_error=0.0;
    *avg_error=0.0;
   }

}  /* brute_reset_distance */
